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Abstract 

LADAR  (LAser  Detection  and  Ranging)  systems  can  be  used  to  provide  2-D 
and  3-D  images  of  scenes.  Generally,  2-D  images  possess  superior  spatial  resolution 
due  to  the  density  of  their  focal  plane  arrays,  but  without  range  data.  A  3-D  LADAR 
system  can  produce  range  to  target  data  at  each  pixel,  but  lacks  the  2-D  system’s 
superior  spatial  resolution.  The  3-D  system  is  limited  by  its  hardware,  specifically  its 
imaging  array.  Gurrently  developers  are  investigating  ways  to  change  the  pixel  size 
in  the  3-D  LADAR  imaging  array,  but  the  cost  of  this  research  is  quite  expensive  and 
technically  challenging.  It  is  the  goal  of  this  work  to  develop  an  algorithm  using  an 
Expectation  Maximization  approach  to  estimate  both  range  and  the  bias  associated 
with  a  3-D  LADAR  system.  The  algorithm  developed  demonstrates  both  spatial  and 
range  resolution  improvement  over  standard  interpolation  techniques  using  both  real 
and  simulated  3-D  and  2-D  LADAR  data. 
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A  Statistical  Approach  to  Fusinc  2-D  and  3-D  LADAR 


Systems 

I.  Introduction 

This  chapter  describes  the  problem  to  be  addressed  by  this  research.  The  back¬ 
ground  of  the  problem  and  goals  for  this  research  are  given,  as  well  as  assump¬ 
tions  used  to  limit  the  scope  of  the  research.  A  discussion  of  previous  related  research 
is  provided  as  well  as  the  organization  for  the  rest  of  the  thesis. 

1.1  Background 

FLASH  3-D  LAser  Detection  and  Ranging  (LADAR)  systems  represent  an  im¬ 
portant  advancement  in  imaging  technology  in  that  they  capture  an  entire  scene 
simultaneously  as  opposed  to  the  way  scanning  systems  form  imagery.  3-D  FLASH 
systems  suffer  from  spatial  resolution  problems  due  to  pixel  pitch  fabrication  limita¬ 
tions.  A  2-D  system  can  produce  high  spatial  resolution  images  but  without  range 
data.  The  3-D  system  can  produce  range  to  target  data,  but  lacks  the  2-D  system’s 
superior  spatial  resolution.  The  3-D  system  is  limited  by  its  hardware,  specifically  its 
imaging  array’s  pixel  pitch.  Pixel  pitch  is  the  distance  between  each  pixel  in  an  array. 
The  pixel  pitch  of  many  3-D  systems  is  100  micro-meters  while  the  2-D  system  pos¬ 
sesses  a  pixel  pitch  of  10  or  25  micro-meters.  Currently  developers  are  investigating 
ways  to  improve  the  pixel  size  in  the  3-D  LADAR  imaging  array,  but  the  costs  of  this 
research  is  quite  expensive  and  technically  challenging. 

Obtaining  better  spatial  resolution  from  a  3-D  LADAR  system  would  improve 
LADAR  systems  with  potential  United  States  Air  Force  (USAF)  LADAR  applica¬ 
tions.  Some  defensive  systems  use  LADAR  for  navigation,  target  recognition,  and 
reconnaissance.  Currently  researchers  are  trying  to  develop  LADAR  systems  that 
help  autonomous  vehicles  navigate  around  unfamiliar  terrain.  One  technology  allows 
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a  3-D  LADAR  system  to  communicate  with  a  vehicles  inertial  measurement  unit  to 
provide  motion  estimates  to  the  vehicle  using  absolute  orientation  [6].  Collections  of 
these  motion  estimates  provide  navigation  to  the  vehicle  as  well  as  other  systems. 

Automatic  target  recognition  is  a  big  focus  for  the  USAF.  This  type  of  technol¬ 
ogy  opens  up  the  use  of  loitering  munitions,  munitions’  miniaturization,  and  target 
recognition.  In  2007  Lockheed  Martin  developed  a  LADAR  based  seeker  head  for  tar¬ 
get  recognition,  known  as  E-LADAR  [13].  This  technology  improves  the  accuracy  of 
many  munitions,  which  improves  the  lethality  of  force  and  decreases  collateral  dam¬ 
age  of  the  munitions.  Pilots  also  use  target  recognition,  a  3-D  target  image  allows 
them  to  make  better  decisions  to  discern  a  target  from  a  non-target.  Improvement  of 
the  system’s  range  estimation  would  make  it  more  accurate  and  would  increase  the 
munitions  resolntion  of  a  target. 

Intelligence  plays  a  big  role  in  today’s  battleheld;  knowing  where  to  go  and 
the  location  of  the  enemy  are  keys  to  maintaining  dominance  of  a  battleheld.  Re¬ 
connaissance  by  LADAR  can  create  3-D  maps  of  whole  scenes  on  the  gronnd  from 
airframes.  After  the  aerial  LADAR  data  is  collected  terrain  maps  are  developed  which 
aid  change  detection  in  battleheld  environments.  Commanders  use  the  3-D  terrain 
maps  in  conjnnction  with  other  intelligence  to  develop  battle  plans  or  investigate  ad¬ 
versaries.  Improving  both  range  estimation  and  spatial  resolntion  of  a  3-D  LADAR 
system  would  allow  autonomous  vehicles  to  hy  closer  to  targets  and  throngh  nrban 
areas,  create  pinpoint  accuracy  for  munitions,  and  improve  battleheld  intelligence  as 
well  as  air  and  space  dominance. 

One  method  for  obtaining  better  spatial  and  range  resolntion  from  3-D  LADAR 
systems  is  to  interpolate  the  images  throngh  various  techniques.  While  this  may 
suffice  as  a  solution  to  the  3-D  pixel  size  problem,  it  is  not  as  accnrate  in  reference 
to  spatial  and  range  resolution  of  the  3-D  images.  Interpolation  may  introdnce  errors 
due  to  aliasing  ehects. 
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1.2  Research  Goals 


The  primary  goal  of  this  research  is  to  prove  that  fusing  both  2-D  and  3-D 
LADAR  images  through  Expectation  Maximization  (EM)  will  increase  the  spatial 
and  range  accuracy  of  the  3-D  system. 

1.3  Assumptions 

For  this  research,  the  following  assumptions  were  made: 

•  The  LADAR  pulse  returns  exist  within  the  range  gate  of  the  system 

•  The  total  LADAR  system  point  spread  function  (PSF)  stays  hxed  and  is  known 

or  can  be  measured 

•  The  LADAR  location  is  known  in  simulated  data 

•  The  2-D  and  3-D  LADAR  system’s  images  are  pre-registered  and  aligned 

1.4.  Related  Research 

This  section  describes  other  methods  to  improve  3-D  LADAR  resolution  and 
range  estimation.  The  methods  discussed  are  interpolation,  microscanning,  texture 
mapping,  and  blind  deconvolution. 

1.4-1  Interpolation.  Interpolation  is  a  method  in  which  new  data  points 
are  created  from  a  set  of  sampled  data  points.  Data  interpolation  can  be  done  in 
numerous  ways,  for  the  purpose  of  this  work  we  will  focus  on  the  pixel  replication 
(zero-order),  linear  (hrst-order) ,  and  cubic  interpolators.  An  interpolator  takes  the 
data  given  to  it  and  creates  new  data  to  a  desired  range.  The  new  data  created  is  an 
approximation  based  on  surrounding  information.  Through  image  interpolation,  each 
pixel  and  its  surrounding  pixels  are  used  to  determine  the  new  pixels.  This  smoothes 
the  data  out  and  allows  one  to  expand  the  resolution  of  an  image.  Interpolating  the 
data  is  a  quick  process  and  is  a  standalone  method.  The  interpolators  only  rely  on 
the  3-D  LADAR  data  and  they  do  not  take  into  account  the  2-D  LADAR  data.  These 
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attributes  make  data  interpolation  an  attractive  method  for  extracting  information 
from  the  3-D  LADAR  system.  Interpolation  does  not  reduce  aliasing  that  already 
exists  in  an  image  [8]. 

1.4- 2  Microscanning.  Microscanning  is  a  method  that  involves  registering 
numerous  3-D  LADAR  data  cubes  in  every  dimension.  It  has  been  proven  to  increase 
spatial  resolution  and  reduce  range  estimation  error  of  3-D  LADAR  systems.  The 
microscanning  process  registers  the  images  of  all  the  data  cubes  at  sub-pixel  resolution 
and  then  uses  interpolation  to  estimate  the  range  of  the  target  [1].  While  under  some 
circumstances  microscanning  may  suffice,  the  system  can  have  some  latency,  due  to 
waiting  on  the  many  frames  it  requires.  If  the  target  ends  up  changing  during  the 
microscanning  process,  frames  could  get  misregistered  and  ultimately  produce  an 
unwanted  image.  The  algorithm  proposed  takes  a  single  3-D  data  cube  and  improves 
its  spatial  resolution  and  range  accuracy. 

1.4.3  Texture  Mapping.  Texture  mapping  is  the  process  of  taking  high 
resolution  2-D  images  and  overlaying  them  onto  3-D  LADAR  scenes  or  aerial  stereo 
imaging.  This  method  improves  the  way  3-D  images  aesthetically  look  but  does  not 
improve  the  range  accuracy  [5].  The  method  registers  the  2-D  image  with  a  3-D 
scene  and  wraps  the  2-D  image  around  the  3-D  scene.  The  data  fusion  proposed  in 
this  research  is  very  different  from  texture  mapping.  The  data  fusion  proposed  uses 
statistics  to  estimate  ranges  based  on  both  3-D  and  2-D  images. 

1.4- 4  Blind  Deconvolution.  There  are  two  previous  research  efforts  that 
used  blind  deconvolution  to  improve  3-D  LADAR  range  estimation,  McMahon’s  [7] 
and  Cain’s  [2]  both  use  a  blind  deconvolution  method.  Cain’s  blind  deconvolution 
method  was  developed  strictly  using  the  Richardson-Lucy  algorithm  to  estimate  the 
pulse  shape  and  then  extract  the  range  from  the  shape.  Cain’s  work  did  not  estimate 
the  intensity  of  the  signals  and  assumed  it  was  known.  McMahon’s  work  is  similar, 
it  uses  the  EM  process  to  come  up  with  an  algorithm  to  estimate  pulse  shape  and 
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extract  the  range  from  it.  McMahon’s  work  deals  with  critically  sampled  data.  The 
data  used  in  this  study  has  been  aliased  to  properly  set  up  a  case  for  3-D  LADAR 
systems,  e.g.  munitions’  seekers,  machine  vision,  etc.  The  proposed  algorithm  also 
estimates  the  range  directly  without  extracting  it  from  the  pulse  shape. 

1 . 5  Thesis  Organization 

Chapter  II  provides  a  description  of  the  LADAR  sensor  model  and  data  de¬ 
velopment  used  for  this  research.  Chapter  III  explains  the  mathematical  derivation 
of  the  algorithm  and  methodology  of  the  research.  Chapter  IV  details  the  results 
from  the  simulations  described  in  Chapter  III.  Finally,  Chapter  V  gives  a  summary  of 
the  research  and  lists  conclusions  of  the  thesis  as  well  as  potential  follow-on  research 
areas. 
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II.  LADAR  Sensor  Model 


This  chapter  provides  the  technical  background  necessary  for  understanding  the 
overall  concepts  of  this  research.  A  description  of  the  LADAR  sensor  model  is 
provided  and  describes  how  data  was  developed  for  this  work. 

2.1  LADAR  Sensor  Models 

2-D  and  3-D  LADAR  systems  interrogate  scenes  through  optics  as  well  as  the 
atmosphere  using  laser  pulses.  Models  exist  for  each  system  that  relate  the  target 
plane  coordinates  {x,y),  to  the  system’s  focal  array  plane  coordinates  {u,v).  The 
coordinates  represent  pixels  in  the  arrays.  The  following  subsections  will  describe  the 
sensor  models  and  how  the  data  was  developed  for  this  work. 

2.1.1  LADAR  Hardware.  Before  discussing  the  details  of  modeling  a  FLASH 
LADAR  system  one  must  understand  its  hardware.  Most  3-D  FLASH  LADAR  sys¬ 
tems  contain  the  same  components  as  shown  in  figure  2.1.  When  a  scene  is  shot,  the 
laser  transmits  through  a  diffuser  (beam  spreader)  in  order  to  cover  the  area  of  the 
scene.  The  lasers  clock  rate  determines  the  pulse  width  of  the  laser.  In  this  research, 
the  pulse  width  of  the  measured  data  is  2.5  nanoseconds.  Using  the  diffuser  ensures 
uniform  illumination  of  the  target.  Once  the  beam  hits  the  target,  the  beam  reflects 
off  the  target  creating  a  return  pulse  back  to  the  system.  This  pulse  is  put  through 
lensing  that  sizes  the  pulse  down  to  the  focal  array  (array  of  receivers).  Focal  array 
sizes  vary  from  system  to  system.  The  light  hitting  the  focal  array  goes  through  a 
data  processor  that  determines  the  ranges  contained  in  the  scene  shot.  The  processed 
data  then  becomes  the  3-D  scene  of  the  target. 

2.1.2  LADAR  sensor  models.  A  2-D  LADAR  system  produces  2-D  intensity 
image  of  a  target.  Most  2-D  systems  can  achieve  a  25  micrometer  or  smaller  pixel 
pitch  which  can  produce  high  resolution  intensity  images.  The  intensity  received 
from  the  target  is  not  the  same  as  the  intensity  of  the  target.  A  relationship  exists 
between  target  intensity  and  the  received  intensity  in  the  focal  plane  and  is  shown  in 
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Figure  2.1:  The  components  of  a  3-D  LADAR  system  while  receiving  a  pulse. 

Eq.  (2.1)  [2].  The  target  intensity  is  A[x,y),  and  the  focal  plane  intensity  is  i{u,v). 
N  is  the  number  of  pixels  in  the  high-resolution  image  plane  in  each  dimension. 


N  N 

i{u,v)  =  ^^A{x,y)h{u-x,v-y)  (2.1) 

x=l  y=l 

Equation  (2.1)  is  a  convolution  of  the  target  intensity  with  the  PSF,  h{x,y). 
This  intensity  model  is  used  to  hnd  the  intensity  of  an  incoherent  system  [4].  The 
LADAR  laser  light  is  coherent  in  nature,  but  its  detected  light  return  is  modeled  as 
incoherent  therefore  the  use  of  Eq.  (2.1)  is  justihed. 

The  PSF  represents  the  effects  of  all  the  optics  as  well  as  the  atmosphere.  The 
optical  transfer  function  (OTF)  of  the  optics,  Hopt{fx,  fy),  is  the  effect  the  optics  have 
on  the  pulse  or  lightwave  coming  through  the  system.  Hoptifx,  fy)  can  easily  found 
by  autocorrelating  the  pupil  function  (t/ens)  as  shown  in  Eq.  (2.3)  [4].  The  variables 
ifx,  fy)  represent  the  spatial  frequencies  in  two  dimensions. 
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^lens{p^)  I/) 


l,for  {x,y)  in  the  aperture 
O,for  {x,y)  outside  the  aperture 


(2.2) 


HoAf^Jy)  = 


X)x=l  X)i;=l  ^le 


,{X  +  y  +  A^)tlensix 


^  flfx  ^  . 

2  ’  y 


A  flfy 


X)x=l  X)i;=l  I'^lensi 


[X, 


(2.3) 


The  variable  A  represents  the  laser  wavelength  and  fi  is  the  focal  length  of  the 
system.  J  represents  the  number  of  pixels  in  each  dimension  of  a  system’s  pupil. 
Hopt{fx,fy)  xn&y  have  range  dependency  but  for  this  research  it  will  stay  fixed  for 
simplicity.  The  atmospheric  OTF,  Hatm{fx,fy)j  is  an  average  and  is  found  using 
Eq.  (2.4)  [10].  The  equation  represents  the  average  effects  of  the  atmosphere  through 
short  exposure.  The  equation  uses  the  wavelength  of  the  laser,  the  focal  length  of 
the  system,  Fried’s  seeing  parameter  (r^),  and  diameter  of  the  lens  {Dr).  Fried’s 
seeing  parameter  is  formed  by  weak  and  strong  atmospheric  turbulence.  A  high 
will  dictate  weak  turbulence,  while  strong  turbulence  is  represented  by  a  small  Tq. 


Hatm{f X:  fy)  6 


-3.44 


1- 


rA2/2(/2+j.2) 


(2.4) 


There  is  also  spatial  blurring  that  is  caused  by  each  detector  due  to  most  im¬ 
ages  being  bigger  than  a  single  detector.  This  blurring  has  its  own  transfer  function 
Hdet{fxi  fy)-  Hdetifx,  fy)  E  Created  from  taking  the  Fourier  transform  of  a  2-by-2  rect 
function  which  is  a  sine  function.  Multiplying  the  optics’  OTF,  atmospheric  OTF,  and 
the  detector  transfer  function  creates  Htot{fxify)  (Eq.  (2.5))  in  the  Fourier  domain. 
Taking  the  inverse  Fourier  transform  of  Htot{fx,fy)  results  in  the  PSF,  h{x,y). 


Htotifx,  fy)  =  Hdetifx,  fy)HatAfx,  fy)Hopt{fx,  fy) 


(2.5) 


The  pulses  from  the  3-D  LADAR  system  will  be  modeled  as  Gaussian  in  time 
as  shown  in  Eq.  (2.6).  By  keeping  A{x,y)  hxed  the  3-D  pulse  becomes  a  function  of 
the  range  at  each  pixel,  r{x,y).  The  variable  represents  discrete  ranges  each  time 
the  pulse  is  sampled  by  the  3-D  system,  for  each  sample  a  2-D  image  is  made  which 
combine  to  make  a  3-D  image.  The  variable  a  is  the  pulse  width  in  meters. 


P{x,y,rk) 


A(x,y)  (r-fc-r(3:,^))^ 

— ) _  e  2^^^ 


(2.6) 


The  3-D  data’s  focal  plane  intensity  is  fonnd  by  convolving  the  PSF,  h{x,y), 
and  the  pulse  as  shown  in  Eq.  (2.7).  Eqnation  (2.7)  also  contains  the  bias,  B{u,v), 
that  is  present  in  the  measured  data.  The  bias  is  generated  is  a  result  of  dark  current 
added  to  the  signal.  This  research  assumes  that  the  bias  follows  a  Poisson  distribution 
due  to  its  discrete  nature  [7].  L  is  the  under  sampling  factor  between  the  2-D  and 
3-D  data.  The  pixel  pitch  of  each  system  determines  the  nndersampling  factor. 


N  N 

I{u,v,rk)  =  EE  P{x,  y,  rk)h{Lu  —  x^  Lv  —  y)  +  B{u,  v)  (2.7) 

x=l  y=l 

The  3-D  data  {d{u,  v,  r^))  is  a  Poisson  random  variable  with  a  mean  of  I{u,  v,  r^). 
The  3-D  data  model  was  chosen  because  incoherent  light  can  be  modeled  as  a  Poisson 
random  variable.  The  LADAR  laser  light  is  coherent  which  makes  the  modeling  of 
the  system  robust.  The  choice  to  treat  the  laser  light  as  incoherent  light  is  based 
on  the  success  of  previous  works  modeling  the  LADAR  light  as  incoherent  [7],  [2]. 
Assnming  every  pixel  and  every  time  instance  in  the  3-D  data  is  independent  the 
joint  probability  {p{d))  of  the  data  is  shown  in  Eq.  (2.8).  The  variable  M  is  number 
of  pixels  in  the  3-D  image  plane  (low  resolution  plane)  in  each  dimension  and  has 
the  relationship  M  =  ^.  The  variable  K  represents  the  total  number  of  2-D  images 
contained  in  the  3-D  data.  Snmming  the  all  the  images  contained  in  the  3-D  data 
creates  the  2-D  images  that  will  be  used  for  the  data  fusion.  The  subsections  below 
discuss  the  two  types  of  data  used  in  this  work. 
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(2.8) 


^  ^  ^  J(u  y  J.  )d2{u,v,rk)^-I{u,v,rk) 

pw=nnn— 


u=l  v=l  k=l 


d2{u,v,rk)\ 


2.1.3  Simulated  data.  The  simulated  images  generated  for  use  in  testing 
the  proposed  algorithm  will  be  50-by-50  (high  resolution)  for  the  2-D  case  and  13- 
by-13  (low  resolution)  for  the  3-D  case.  These  resolutions  produce  an  undersampling 
factor  of  approximately  four  (L  4).  The  undersampling  factor  and  resolutions 
were  chosen  because  currently  2-D  LADAR  systems  can  easily  achieve  a  pixel  pitch 
of  25  micrometers  while  3-D  LADAR  systems  possess  a  100  micrometer  pixel  pitch 
which  produces  an  undersampling  factor  of  four.  The  simulated  3-D  data  contains  20 
simulated  pulse  returns. 

The  hrst  target  used  to  produce  the  simulated  data  is  shown  in  hgure  2(a). 
The  hgure  depicts  a  two  building  target  that  is  imaged  onto  a  50-by-50  target  plane. 
The  top  of  the  buildings  are  located  at  10000  meters  and  the  ground  is  located  at 
10002  meters.  The  raw  under  sampled  3-D  data  produced  an  estimated  range  shown 
in  hgure  2(b).  The  estimated  range  is  the  result  of  using  a  cross-correlation  range 
estimator  (matched  hlter)  which  will  be  discussed  in  detail  in  section  3.2. 


Figure  2.2:  (a)  The  50-by-50  target  area,  (b)  The  13-by-13  raw  estimated  range 

from  the  undersampled  3-D  data  on  the  low  resolution  grid. 
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The  3-D  data  can  also  be  represented  by  its  2-D  images  at  each  time  instance. 
A  time  instance  is  governed  by  the  pulse  width  of  the  laser.  For  the  simulated  systems 
the  pulse  width  is  two  nanoseconds.  The  pulse  width  dictates  the  time  it  takes  the 
pulse  to  hit  the  target  and  return  back  to  the  receiver.  Figure  2.3  shows  the  2-D 
images  of  the  data  cube  for  every  other  time  instance.  The  images  represent  the 
image  received  by  the  LADAR  system  every  four  nanoseconds. 


Figure  2.3:  2-D  intensity  images  of  the  3-D  data  cube  over  the  course  of  36  nanosec¬ 

onds.  Each  image  is  shown  on  a  13-by-13  pixel  grid  and  represents  the  data  slice  taken 
every  4  nanoseconds. 


An  enlarged  image  at  the  second  time  instance  is  shown  in  figure  2.4,  the  figure 
includes  both  high  and  low  resolution  images.  The  low  resolution  image  is  what  a  3-D 
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LADAR  system  provides  in  2-D.  The  sampling  in  this  example  makes  it  difficult  to 
discern  the  actual  structure  of  the  object.  The  3-D  data  is  13-by-13  and  is  pictorially 
different  from  the  original  target  area.  This  provides  the  motivation  for  trying  to 
improve  the  3-D  data  through  various  methods.  Figure  2.3(c)  shows  the  2-D  image 
intensity  which  represents  what  a  2-D  camera  would  see.  The  image  was  created  by 
summing  the  high  resolution  images  in  the  high  resolution  3-D  data  cube. 


(c) 


Figure  2.4:  (a)  The  low  resolution  (13-by-13)  data  slice,  (b)  The  high  resolution 

(50-by-50)  data  slice,  (c)  The  2-D  data  to  be  fused  with  the  3-D  data. 

The  second  target  used  to  produce  the  simulated  data  is  more  complicated. 
The  target  contains  numerous  buildings  at  different  heights  while  the  ground  is  still 
located  10002  meters  from  the  LADAR  system.  The  second  target’s  profile  is  shown 
in  figure  2.5. 
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(a)  (b) 


Figure  2.5:  (a)  50-by-50  target  area  of  the  second  simulation,  (b)  The  2-D  image 

resulting  from  summing  the  second  targets  high  resolution  3-D  data  cube. 


Figure  2.6:  Setup  to  obtain  measured  data. 


2.1.4  Measured  data.  The  measured  data  was  obtained  from  a  3-D  LADAR 
system  that  shot  a  scene  of  a  3  bar  cut  out  (third  target)  as  depicted  in  figure  2.6. 
The  system  used  17  time  samples  to  build  a  3-D  data  cube.  The  system  possesses  a 
128-by-128  grid  of  detectors.  The  system  created  20  separate  3-D  data  cubes.  Each 
image  in  the  cube  was  cropped  down  to  64-by-64  in  order  to  focus  on  the  target  area. 

Unlike  the  simulated  data,  the  measured  data  exhibits  non-ideal  behavior  which 
includes  but  is  not  limited  to  gain  variation  and  electronic  noise.  When  the  LADAR 
system  shoots  a  scene  a  gain  drop  exists  due  to  the  laser  being  incident  on  a  larger 
area  in  the  detector  array  [11].  This  phenomena  causes  a  gain  variation  between 
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the  first  and  second  snrfaces.  The  LADAR  system  used  to  create  the  data  is  not  a 
photon  counting  system.  In  photon  counting  systems  every  detected  return  equals  one 
photon.  The  proposed  algorithm  does  not  take  into  account  gain  variation  or  photon 
counting  conversion  factors  in  order  simplify  the  algorithms  development.  [11]  shows 
how  to  correct  for  this  behavior.  The  data  used  is  raw  data  that  contain  these  effects 
in  it  without  correction. 

Another  parameter  is  needed  to  use  the  measured  data  and  it  is  the  system’s 
optical  cutoff  frequency,  fc-  The  cutoff  frequency  of  the  system  is  determined  by  using 
Eq.  (2.9)  [4].  The  diameter  [D)  and  focal  distance  {Z)  are  shown  in  hgure  2.6  and 
the  wavelength  (A)  of  the  system  is  1.55  micrometers.  In  order  to  achieve  the  correct 
sampling  rate  the  Nyquist  rate  must  be  used  which  is  two  times  the  cutoff  frequency. 
This  makes  the  sampling  rate  116.25  micrometers  [10].  Since  The  pixel  pitch  is  100 
micrometers,  the  3-D  cube  from  the  data  is  oversampled. 


fc  = 


D 

\Z 


(2.9) 


Since  the  measured  data  is  critically  sampled  the  true  range  can  be  extracted 
from  it.  The  true  range  is  created  from  averaging  all  the  ranges  from  each  data  cube 
and  subtracting  the  mean  of  the  analyzed  cube.  The  hfth  data  cube  was  used  for  this 
research.  The  measured  data  range  information  is  in  terms  of  pulse  returns  (samples). 
The  true  range  for  the  measured  data  is  shown  in  hgure  2.7. 

Downsampling  the  measured  data  creates  under  sampled  data  that  will  be  used 
for  analysis  in  this  research.  This  downsampling  effect  captures  the  effect  of  a  larger 
pixel  and  data  aliasing.  The  effect  of  downsampling  an  image  in  the  cube  is  shown  in 
hgure  2.8. 

The  downsampling  ehect  shown  in  hgure  2.6(a)  is  also  seen  in  all  the  images  of 
the  measured  data  cube.  The  images  produced  by  every  other  pulse  return  are  shown 
in  hgure  2.9. 
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Figure  2.7:  Target  area  of  measured  data  on  a  64-by-64  grid. 
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(a)  (b) 


t 


(c) 

Figure  2.8:  (a)  The  low  resolution  (16-by-16)  data  slice,  (b)  The  high  resolution 

(64-by-64)  data  slice,  (c)  The  2-D  data  that  will  be  fused  with  the  3-D  data. 
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Figure  2.9:  2-D  intensity  images  of  the  3-D  measured  data  cube.  Each  image  is 

shown  on  a  16-by-16  pixel  grid. 


17 


III.  Research  Methodology 


This  chapter  describes  the  methods  used  to  derive  the  proposed  sensor  fusion  al¬ 
gorithm  as  well  as  interpolation  methods  used  for  comparison.  The  first  section 
will  describe  the  use  of  the  EM  algorithm  to  estimate  both  the  range  and  bias  of  the 
LADAR  signals.  The  second  section  will  detail  the  types  of  interpolation  that  will  be 
used  for  comparison  to  the  proposed  algorithm.  Using  the  equations  in  simulations 
and  the  results  of  those  simulations  will  be  the  subject  of  chapter  4. 


3.1  Proposed  Sensor  Fusion  Algorithm 

The  proposed  algorithm  is  derived  using  an  EM  approach  to  estimate  the  range 
r{x,y)  to  each  pixel  in  the  target  area  and  bias  B{u,v)  for  each  detector  in  the  array. 
The  proposed  sensor  fusion  approach  will  estimate  A{x,  y)  by  using  a  Richardson-Lucy 
deconvolution  [9]  of  the  properly  sampled  2-D  data  as  shown  in  Eq.  (3.1).  This  tech¬ 
nique  is  an  iterative  process  in  which  each  iteration  produces  an  Anew(x,  y).  Aoid(x,  y) 
represents  the  estimate  for  A[x,y),  d2d{u,v)  is  the  observed  2-D  data,  and  i{u,v)  is 
the  estimated  focal  plane  intensity.  h{x,y)  represents  the  total  optical  point  spread 
function  (PSF).  For  this  work,  it  is  assumed  that  the  PSF  is  known.  The  initial 
estimate  for  A{x,y)  is  a  constant. 


A(^X,  y^new 


N  N 

Mx,y)oid^Yl 

v=l  u=l 


d2d{u,v) 
i{u,  v) 


h{u 


x,v-y) 


(3.1) 


The  GEM  approach  proposed  by  McMahon  [7]  is  similar  to  the  EM  approach 
in  this  work.  They  are  not  the  equivalent  because  McMahon’s  work  does  not  deal 
with  undersampled  data  and  does  not  feature  the  use  of  2-D  and  3-D  data.  The  first 
step  of  the  EM  approach  is  to  create  a  statistical  model  for  the  measured  data,  which 
is  known  as  the  incomplete  data.  Inventing  a  set  of  mythical  data  (complete  data) 
and  its  relationship  to  the  incomplete  data  is  the  second  step.  The  third  step  is  to 
select  a  statistical  model  for  the  complete  data  such  that  it  adheres  to  the  relationship 
of  the  complete  to  incomplete  data.  Next  is  to  form  a  complete  data  log-likelihood. 
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In  step  five,  the  conditional  expectation  of  the  complete  log-likelihood  is  computed 
with  respect  to  the  incomplete  data.  The  last  step  is  to  maximize  the  conditional 
expectation  with  respect  to  the  parameter  that  is  being  estimated.  The  following 
subsections  describe  the  EM  approach  to  estimate  the  range  [3]. 


3.1.1  Formulating  complete  and  incomplete  data.  As  stated  in  chapter  2, 
the  3-D  LADAR  data  (incomplete  data)  is  a  realization  of  a  Poisson  random  variable 
at  each  pixel.  The  joint  probability  of  the  incomplete  data  is  shown  in  Eq.  (2.8).  Due 
to  the  bias  estimation  the  complete  data  will  contain  two  variables,  di{u,v,x,y,rk) 
and  d2{u,v,rk).  The  incomplete  data,  d,  has  a  relationship  to  the  complete  data  as 
shown  in  Eq.  (3.2). 

M  M 

d{u,v,rk)  =  EE  di{u,v,x,y,rk)  +  d2{u,v,rk)  (3.2) 

x=l  y=l 


The  variables  {x,  y)  represent  the  target  plane  pixel  locations  and  {u,  v)  repre¬ 
sents  the  focal  plane  coordinates.  The  variable  represents  discrete  ranges  each  time 
the  pulse  is  sampled  by  the  3-D  system.  The  expectations  (means)  of  the  complete 
data  are  shown  in  Eq.  (3.3).  The  means  were  chosen  based  on  the  complete  data’s 
relationship  to  the  incomplete  data. 


E[di{u,v,x,y,rk)] 


E[d2{u,v,rk)] 


B{u,v) 


x,Lv  -  y) 


Since  the  sum  of  two  Poisson  random  variables  is  Poisson,  the  complete  data 
can  be  modeled  as  a  Poisson  random  variable  at  each  pixel,  retaining  the  validity  of 
Eq.  (3.3).  Based  on  the  chosen  means,  the  complete  data  probabilities  pi  and  p2  are 
shown  in  Eq.  (3.3). 
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h{Lu—x,Lv—y) 


Pi[{di{u,v,x,y,rk)]  = 

h{Lu  -X,Lv-  y)dii^X,x,y,rk)^-^, 


.  (rfc-r(a;,y))^ 


di{u,v,x,y,rky. 


P2[d2{u,v,rk)]  = 

B[u,  vY2iu,v,rk) B(u,v) 

d2{u,v,rk)\ 


(3.3) 


Taking  the  product  of  the  probabilities  at  each  sample  location  (x,  y,  u,  v)  and 
each  range  (r^)  produces  the  joint  probability  (ps)  for  the  complete  data,  Eq.  (3.4). 


P3  = 


'  M  M  K  M  M 

nnnnn  pi[{di{u,v,x,y,rk)\ 

u=l v=l k=l x=l y=l 


X 


M  M  K 


nnn  P2[d2{u,v,rk)\ 


u=l  v=l  k=l 


(3.4) 


Taking  the  natural  log  of  Eq.  (3.4)  produces  the  complete  data  log-likelihood 
shown  in  Eq.  (3.5). 


M  M  M  M  K 


i  =  E  E  E  E  E  l*(«.  V,  nM^e'^^HLu  -  X,  Lv  -  y)) 

x=l  y=l  u=l  v=l  k=l  ' 


A{x,y)  (r^-rixy)P 

-e  2^^  h{Lu  —  X,  Lv  —  yj  —  ln[di{u,  v,  x,  p,  rk)'.)] 


\fLna 

M  M  K 

EEE  [d2{u,  V,  rk)ln{B{u,  v))  -  B{u,  v)  -  ln{d2{u,  v,  r^)!)]  (3.5) 


u=l  v=l  k=l 


3.1.2  Finding  the  expectation.  The  conditional  expectation  of  the  complete 
data  log-likelihood  is  found  by  taking  the  expectation  with  respect  to  the  incom¬ 
plete  data,  the  pulse  estimate  {Poid{x,y,rk)),  and  the  bias  estimate  {Boid{u,v)).  The 
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conditional  expectation  (Q)  in  general  form  is  shown  in  Eq.  (3.6). 


M  M  M  M  K 

EEEEE  [E[di{u,  V,  X,  y,  rk)\d{u,  v,  Vk),  Poid{x,  y,  r^),  Boid{u,  v)] 

x=l  y=l  u=l  v=l  k=l 

A(x^  tf]  '^(^1?/))^ 

X  (/n( —  '  h{Lu  —  X,  Lv  —  y))  +  /n(e  2^:^  ))  —  P{x,  y,  rk)h{Lu  —  x,  Lv  —  y) 

y2Tia 

-E[ln{di{u,  V,  X,  y,  rky.)\d{u,  v,  r^),  Poid{x,  y,  rk),  Boid{u,  n)]] 

M  M  K 

+EEE  [E[{d2{u,  V,  rk)\d{u,  n,  rfc),  Poid{x,  y,  rk),Boid{u,  v)]ln{B{u,  v)) 

u=l  v=l  k=l 

-B{u,v)  -  E[{d2{u,v,rky.\d{u,v,rk),Poid{x,y,rk),Boid{u,v)]]  (3.6) 


The  variable  loid  is  the  image  produced  from  the  pulse  estimate,  Eq.  (3.7),  which 
is  needed  to  hnd  each  of  the  individual  conditional  expectations. 

N  N 

Iold{Ui  X,  Tk)  =  EE  Poid{x,  y,  rk)h{Lu  -  x,Lv  -  y)  +  Boid{u,  v)  (3.7) 

x=l  y=l 

The  solution  for  the  individual  conditional  expectations  are  shown  in  Eq.  (3.8) 

[12]. 


E[di(u,  V,  X,  y,  rk)\d(u,  v,  rk),  Poid(x,  y,  Vk),  Boid(u,  n)] 
d(u,  V,  rk)Poid(x,  y,  rk)h(Lu  -x,Lv  -  y) 
Ioid(u,v,rk) 

E[d2(u,  V,  rk)\d(u,  v,  rk),Poid(x,  y,  Vk),  Boid(u,  n)] 


d(u,v,rk)Boid(u,v) 

Ioid(u,v,rk) 


(3.8) 
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Substituting  the  previous  result  back  into  Q  produces  Eq.  (3.9). 


M  M  M  M  K 

ZZEEEt 

x=l  y=l  u=l  v=l  k=l 


X  {ln{ — h{Lu  —  x,Lv  —  y))  — 
v27r(T 

A{x,y)  Ufc-r-qt/))" 

-e  y[Lu  —  x,Lv  —  y) 


d{u,  V,  rk)Poid{x,  y,  rk)h{Lu  -x,Lv  -  y) 
Ioid{u,v,rk) 

irk-r{x,y)f 


2a2 


vra 


-E[ln{di{u,  V,  X,  y,  rfc)!)  |d(M,  v,  Vk),  Poid{x,  y,  Vk),  Boid{u,  u)]] 


M  M  K 


rd{u,v,rk)Boidiu,v) 

-EEE  [  I  ,Ju  V  n) 

u=iv=ik=i  ioid{u,u,rk) 

-B{u,  v)  -  E[{d2{u,  V,  rk)\\d{u,  v,  rk),Poid{x,  y,  r^),  Boid{u,  u)]] 


(3.9) 


3.1.3  Maximizing  the  Expectation.  With  the  conditional  expectation  found 
the  next  step  is  to  maximize  it  with  respect  to  the  range  and  the  bias.  Taking  the 
derivative  of  Q  with  respect  to  r{xo,yo),  then  setting  the  derivative  equal  to  zero 
and  solving  for  r{x,y)  will  maximize  the  conditional  expectation  for  the  range  at 
each  point  in  the  scene.  Doing  the  same  with  respect  to  the  bias  will  maximize  the 
conditional  expectation  for  the  bias.  The  bias  and  range  will  be  estimated  separately. 
In  order  to  estimate  the  range  an  assumption  is  made  that  the  pulse  always  exists 
in  the  range  gate.  This  removes  all  range  dependence  from  the  sum  of  I{u,v,rk)  as 
shown  in  Eq.  (3.10).  The  variable  C  represents  a  constant. 


M  M  M  M  K 


(3.10) 


x=l  y=l  u=l  v=l  k=l 


Tia 


Since  the  bias  portion  (second  summation  group  of  Eq.  (3.9))  of  Q  does  not 
depend  on  r{x,y),  its  derivative  with  respect  to  the  range  goes  to  zero.  Given  this 
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property  and  Eq.  (3.10),  Q  reduces  to  Eq.  (3.11). 


M  M  M  M  K 

ZEEEEi 

x=l  y=l  u=l  v=l  k=l 


d{u,  V,  rk)Poid{x,  y,  rk)h{Lu  -x,Lv  -  y) 
Ioid{u,v,rk) 

(rfc  -  r{x,y)y 


-  X,  Lv  -  y))  -  , 

v27r(T  2(7^ 

-C  -  E[ln{di{u,v,x,y,rky.)\d{u,v,rk),Poid{x,y,rk),Boid{u,v)]]  +  0  (3.11) 


Taking  the  derivative  of  the  hrst  piece  of  Eq.  (3.11)  with  respect  to  r{xo,yo) 
results  in  Eq.  (3.12) 


d  d{u,  V,  rk)Poid{x,  y,  rk)h{Lu  -x,Lv  -  y) 


d{r{xo,yo)) 


Ioid{u,v,rk) 


T  w  {rk-r{x,y)f 

X  (/n(  h{Lu  -x,Lv-  y)) - - )]  = 

y/zna  za 

d{u,v,rk)Poid{x,y,rk)h{Lu  -  x,Lv  -  y)  d  -{rk-r{x,y)Y  _ 

Ioid{u,v,rk)  d{r{xo,yo))  2a‘^ 

d{u,  V,  rk)Poid{x,  y,  rk)h{Lu  -x,Lv-  y)  {2rk  -  2r{x,  y))  d{r{x,  y)) 


Ioid{u,v,rk)  20-2  d{r{xo,yo)) 

d{u,  V,  rk)Poid{x,  y,  rk)h{Lu  -x,Lv-  y)  (r^  -  r{x,  y)) 


Iold{u,  V,  Tk) 


)5{x  -  Xo,y  -  yo)  (3.12) 


The  derivative  of  the  second  piece  of  Eq.  (3.11)  is  shown  to  be  zero  in  Eq.  (3.13). 
This  result  is  due  to  the  conditional  expectation  of  the  complete  data,  given  the 
incomplete  data  and  old  estimates,  will  not  be  a  function  of  the  new  estimates  thus 
rendering  zero  dependency  on  r{x,y). 


d 


d{r{xo,yo)) 


E[ln{di{u,v,x,y,rky.)\d{u,v,rk),Poid{x,y,rk),Boid{u,v)]  =  0  (3.13) 
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Using  Eqs.  (3.10),  (3.12),  and  (3.13),  the  general  form  is  shown 


in 


Eq.  (3.14). 


dQ 


d{r{xo,yo)) 

M  M  M  M  K  r 

ZZZEE 

x=l  y=l  u=l  v=l  k=l 

(rfc  -  r{x,y)) 


d{u,v,rk) 


X  ■ 


a 


Ioidiu,v,rk) 

2  6{x  -  xo,y  -  yo) 


Poid{x,  y,  rk)h{Lu  -x,Lv  -  y) 


(3.14) 


Applying  the  sifting  property  to  Eq.  (3.14)  removes  the  summations  over  x  and 
y  and  results  in  the  Q  that  will  be  maximized  (Eq.  (3.15)). 


dQ 


2^  2^  2^  T  ..  ^  -.PoidAo,  I/O,  rk)h{Lu-xo,  Lv-yo) - 


d{r{xo,yo))  Ioid{u,v,rk) 

Setting  Eq.  (3.15)  equal  to  zero  results  in  Eq.  (3.16). 


(3.15) 


M  M  K 


0  =  E  E  E  2T'^'‘\Po>Ax„,y,,r,)h{Ln  -  x„.  Lv  -  mA"  ~ 

,,j-Aloid{u,v,rk) 

u=l  v=l  k=l  \  ^  ^  / 


Distributing  and  r{xo,yo)  makes  Eq.  (3.16)  easier  to  manipulate  (Eq.  (3.17)). 


0  = 


M  M  K  X 

ixo,yo)  EEEivra  Poid{xo,  yo,  rk)h{Lu  -  xo,  Lv  -  yo) 

u=l  v=l  k=i  Pid[u,v,rk) 


M  M  K 


d{u,v,rk) 


-  y]  y]  y]  rk^^^^^Y~^^^Poid{xo,  yo,  rk)h{Lu  -  Xo,  Lv  -  yo)  (3.17) 

Moving  the  terms  from  one  side  to  the  other  in  Eq.  (3.17),  then  separating  the 
terms  provides  a  solution,  Eq.  (3.18),  that  is  iterative  and  updates  r{xo,yo)  for  each 
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iteration. 


_  Ef=i  rkPoid{xo,  yo,  Tk)  Eili  Ejli  tSSS)  “  ^0,  Lv  -  yo) 

^{Xo^yQ)  —  ^  ^  ^  ^  ,-^M  ■s-^M  d('ii..i,.r,A  ,/^  ^  X  t'J-toJ 


ELl  Poldixo,  yo,  Tk)  E!1i  E!!1i  tSSS)  “  ^0’ 


The  same  process  is  used  to  solve  for  the  bias.  The  range  portion  (first  summa¬ 
tion  group  of  Eq.  (3.9))  of  Q  does  not  contain  the  bias  variable  and  does  not  depend 
on  Q,  its  derivative  with  respect  to  B{uo,vo)  is  zero  reducing  Q  to  Eq.  (3.19). 


Q  = 


M  M  K 


rd{u,v,rk)Boidiu,v) 

Z.  2.  Z.  [  I  ,Ju  V  n) 

„=i  t,=i  fc=i  ioid{u,u,rk) 

-B{u,  v)  -  E[{d2{u,  V,  rk)\\d{u,  v,  rk),Poid{x,  y,  r^),  Boid{u,  n)]] 


(3.19) 


Taking  the  partial  derivative  the  hrst  piece  of  Eq.  (3.19)  with  respect  to  the 
bias,  B{u,v),  is  shown  in  Eq.  (3.20). 


d  d{u,v,rk)Boid{u,v) 
d{B{uo,vo))  Ioid{u,v,rk) 
d{u,v,rk)Boid{u,v)  d 


ln{B{u,v))  = 
ln{B{u,v))  = 


Ioidiu,v,rk)  d{B{uo,Vo)) 

d{u,v,rk)Boid{u,v)  d{B{u,v)) 
Ioid{u,  V,  rk)B{u,  v)  d{B{uo,  Vq)) 
Boid{u,v)d{u,v,rk) 


Ioid{u,v,rk)B{u,v 


-6{u-Uo,v  -  Vo) 


(3.20) 


The  partial  derivative  of  B{u,v),  second  piece  of  Eq.  (3.19),  results  in  a  Dirac 


delta. 


d 


d{B{uo,vo)) 


B{u,  v)  =  S{u  —  uo,v  —  no) 


(3.21) 
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The  derivative  of  the  last  piece  of  Eq.  (3.19)  is  shown  to  be  zero  in  Eq.  (3.22). 
This  result  is  due  to  the  conditional  expectation  not  depending  on  new  estimates. 


d 


d{B{uo,vo)) 


E[ln{d2{u,v,rky.)\d{u,v,rk),Poid{x,y,rk),Boid{u,v)]  =0  (3.22) 


Using  Eqs.  (3.20),  (3.21),  and  (3.22),  the  general  form  of  the  derivative  of  Q 
with  respect  to  B{uq,Vq)  becomes  Eq.  (3.23). 


dQ  ^ 
d{B{uo,vo)) 

M  M  K 

-  uo,v  -  Vo) 


Boid{u,v)d{u,v,rk) 

Ioid{u,v,rk)B{u,v) 


1 


(3.23) 


Applying  the  sifting  property  to  Eq.  (3.23)  results  in  Eq.  (3.24) 
dQ  _  Boid{uo,vo)d{u,v,rk)  _  ^ 

d{B{uo,  Vo))  ^  Ioid{u,  V,  rk)B{uo,  vo) 


(3.24) 


Setting  Eq.  (3.23)  equal  to  zero  and  separating  terms  results  in  the  direct  solu¬ 
tion  for  estimating  the  bias,  Eq.  (3.25).  The  solution  is  iterative  and  updates  B{uo,  vo) 
for  each  iteration.  The  solutions  for  both  range  and  bias  are  iterative  and  must  be 
solved  concurrently  to  estimate  each  one  properly. 

B{uo,  Vo)  =  Bm{uo,  Vo)  ^3  25) 


From  the  equations  described  in  this  section,  the  algorithm  begins  by  estimating 
properly  sampled  2-D  data.  The  algorithm  then  uses  the  2-D  and  3-D  data  to  estimate 
both  range  and  bias.  The  estimates  are  then  used  to  create  new  estimates  with 
each  iteration.  This  provides  for  improved  range  and  bias  estimates  every  iteration, 
ultimately  improving  the  range  accuracy  of  the  raw  3-D  data. 
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3.2  Data  Interpolation 

Interpolation  is  a  method  in  which  new  data  points  are  created  from  a  set  of 
under  sampled  data  points.  Data  interpolation  can  be  done  in  numerous  ways,  for 
the  purpose  of  this  paper  we  will  focus  on  the  pixel  replication  (zero-order),  linear 
(hrst-order) ,  and  cubic  interpolators  [8].  An  interpolator  takes  the  data  given  to  it 
and  creates  new  data  to  a  desired  range.  In  this  paper,  the  interpolator  will  take  a 
13-by-13  image  and  create  a  50-by-50  image.  Interpolating  the  data  is  a  relatively 
quick  process  and  is  a  standalone  method.  These  attributes  make  data  interpolation 
an  attractive  method  for  extracting  information  from  the  3-D  LADAR  system. 

Each  interpolator  has  a  single  basis  equation,  which  is  manipulated  at  each 
higher  order.  The  basis  equation  is  representative  of  what  an  interpolator  performs 
on  image  i{n,m),  an  A^i-by-Mi  under  sampled  image.  Eq.  (3.26)  shows  the  basis 
equation  for  interpolating  data.  The  first  part  of  the  equation  is  creating  a  comb 
of  the  image,  icomb,  that  is  sampled  at  rate  L.  L  is  determined  by  the  following 
equality,  MiL  =  M.  For  example,  producing  a  50-by-50  image  from  a  25-by-25  image 
would  dictate  L  =  2.  The  comb  spreads  the  under  sampled  image  out  to  the  desired 
resolution,  placing  zeros  around  each  pixel.  A  sample  comb  and  image  pair  is  shown 
in  hgure  3.1.  Figure  1(b)  shows  how  the  comb  spreads  the  under  sampled  image  out 
to  the  desired  resolution,  placing  zeros  around  each  pixel.  Convoluting  the  comb  with 
a  hlter,  g{n,  m),  produces  the  interpolated  image.  The  hlter  type  dictates  the  type  of 
interpolator  being  used. 


Ni  Ml 

ic(ymb{P'2i^2)  EE  i(n,  m)5{Ln  —  n2)S{Lm  —  m2) 

n=l  m=l 

N  M 

iint{n,m)  =  icomb{n2,m2)g{n  -  n2,m  -  m2)  (3.26) 

n2=l  m2=l 

For  the  case  of  the  pixel  replication  interpolator  the  hlter  used  is  a  rect  that 
is  L-hj-L.  This  means  that  the  pixel  replication  interpolator  takes  the  closest  pixels 
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O  lA  Q 


Figure  3.1:  (a)  32-by-32  sample  image,  (b)  Comb  of  the  sample  image  on  a  128 

pixel  grid. 

and  replicates  them  out  to  L  pixels  around  that  pixel.  The  linear  interpolator  uses 
a  triangle  hlter  that  is  created  by  convoluting  two  L-hj-L  rect  hlters.  This  allows 
for  the  use  of  another  pixel’s  information  in  the  area  of  the  interpolated  pixel,  rather 
than  just  replicating  it.  Finally,  the  cubic  interpolator  uses  a  piecewise  function  as 
its  hlter.  Eq.  (3.27)  shows  the  piecewise  function  in  the  x-dimension.  The  coefficient 
a  is  calculated  based  on  the  interpolation  size,  N-hj-M. 


{a  +  2)  |a;|^  —  (a  +  3)  \x\^  +  1,  0  <  |a;|  <  1 


9i.x)  =  { 


a  \x\ 


5a  \x\ 


4a, 


1  <  b|  <  2 


0,  ‘2  <  \x 


(3.27) 


The  interpolators  only  rely  on  the  3-D  LADAR  data  and  they  do  not  take  into 
account  the  2-D  LADAR  data.  Interpolation  is  not  complete  itself  without  a  way 


to  estimate  the  range.  This  research  will  apply  a  cross-correlation  range  estimator 
(matched  hlter)  to  all  the  interpolated  images.  The  matched  hlter  uses  a  cross  corre¬ 
lation  function  shown  in  Eq.  (3.28)  [10]  and  assumes  a  target  is  detected.  The  C{R) 
represents  the  correlation  of  the  range,  Dk  is  the  data,  Pt[k  —  2R/ (cA  t)]  is  the  wave 
form.  The  value  R  represents  a  single  range  in  a  set  of  ranges.  The  range  set  deter¬ 
mines  how  hne  the  range  estimates  become.  The  function  produces  a  value  for  each 
range  in  the  set,  the  range  affiliated  with  the  maximum  value  is  the  range  chosen. 

C{R)  =  DkPt[k  -  2R/{cA  t)]  (3.28) 

k=l 
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IV.  Results  and  Analysis 


This  chapter  presents  the  results  of  applying  the  methods  described  in  Chapter  3 
to  the  data  described  in  Chapter  2.  Section  4.1  details  the  results  found  when 
all  four  methods  are  applied  to  the  simulated  data.  Section  4.2  discusses  the  results 
of  the  methods  applied  to  the  measured  data.  The  root  mean  square  error  (RMSE) 
and  graphs  of  the  range  will  represent  the  results  from  each  method.  This  analysis 
will  calculate  the  RMSE  between  the  target  area  (true  range)  and  estimated  ranges 
as  shown  Eq.  (4.1).  Rangcest  is  the  estimated  range  of  the  reconstruction  method  and 
Rangetruth  is  the  target  area  described  in  Chapter  2.  The  variable  M  represents  the 
number  of  pixels  in  each  dimension  on  the  high  resolution  grid. 


RMSE 


Eili  Y.^=iiR(^^9eest{x,y)  -  Rangetruth{x,y)Y 

M2 


(4.1) 


4.-1  Simulated  target  results 

The  first  target  discussed  in  Chapter  2  was  put  through  400  iterations  of  the  pro¬ 
posed  algorithm,  starting  at  a  flat  range.  After  400  iterations,  the  algorithm  achieved 
a  root  mean  square  error  of  approximately  222  milimeters.  Figure  4.1  demonstrates 
how  fast  the  update  moved  to  the  achieved  range  RMSE.  The  number  of  iterations 
where  chosen  based  on  the  graph  shown  in  figure  4.1,  the  full  convergence  of  the 
algorithm  happens  at  approximately  iteration  325. 

After  estimating  the  range  of  the  first  target  with  the  proposed  algorithm,  the 
Pixel  Replication  interpolator  was  used  on  the  data.  The  Pixel  Replication  interpo¬ 
lator  takes  the  data  and  spreads  it  out  from  13-by-13  grid  of  pixels  to  50-by-50  grid. 
This  interpolator  achieved  a  range  RMSE  of  333  milimeters.  Then  the  linear  interpo¬ 
lator  was  applied  to  the  data  as  well.  The  linear  interpolator  achieved  a  range  RMSE 
of  427  milimeters.  The  cubic  interpolator  achieved  a  range  RMSE  of  420  milime¬ 
ters.  The  estimated  ranges  produced  by  the  interpolators  and  proposed  algorithm  are 
shown  in  figure  4.2.  Figure  4.2  shows  the  range  portrayed  on  to  a  3-D  grid.  Table 
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Figure  4.1:  Range  RMSE  vs  iterations  after  applying  the  proposed  algorithm  to 

the  hrst  target.  This  graph  shows  the  RMSE  reducing  drastically  between  zero  and 
hfty  iterations.  Due  to  the  steady  tail  of  the  data,  more  than  400  iterations  of  the 
algorithm  would  only  make  minor  reductions  in  RMSE  of  the  range. 

4.1  shows  the  comparison  of  all  methods  of  improving  range  estimation  for  the  hrst 
target. 

The  proposed  algorithm  then  used  to  estimate  the  range  to  the  second  target 
starting  from  a  hat  range.  Due  to  the  beam  shape  the  error  was  taken  within  a  36- 
by-28  box  starting  from  pixel  location  (10, 14).  The  number  of  iterations  completed 
by  the  proposed  algorithm  was  800  and  resulted  in  a  RMSE  of  132  milimeters.  Figure 
4.3  shows  the  RMSE  as  a  function  of  iteration  count. 

Interpolators  were  then  used  to  estimate  the  range  of  the  second  target.  The 
pixel  replication  interpolator  achieved  a  range  RMSE  of  417  milimeters.  The  lin¬ 
ear  interpolator  achieved  a  range  RMSE  of  451  milimeters.  The  cubic  interpolator 


31 


Figure  4.2:  (a)  50-by-50  range  truth  of  the  hrst  target.  First  target  50-by-50 

estimated  ranges:  (b)  Pixel  replication  estimated  range,  (c)  Linear  estimated  range, 
(d)  Cubic  estimated  range,  (e)  EM  algorithm’s  estimated  range. 
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Table  4.1:  Method  comparison  for  the  first  target 


Method 

Range  RMSE 

Proposed  Algorithm 
Pixel  Replication 
Linear  Interpolation 
Cubic  Interpolation 

222  mm 

333  mm 

427  mm 

420  mm 

achieved  a  range  RMSE  of  445  milimeters.  Due  to  the  high  error  of  each  interpolator’s 
estimated  range,  the  visual  result  of  the  range  estimates  will  be  shown  on  a  50-by-50 
2-D  range  map  and  are  shown  in  hgure  4.4.  The  color  bar  in  each  subhgure  represents 
range  in  meters.  Table  4.2  shows  the  comparison  of  all  methods  of  improving  range 
estimation  for  the  second  target. 


Table  4.2:  Method  comparison  for  the  second  target 


Method 

Range  RMSE 

Proposed  Algorithm 
Pixel  Replication 
Linear  Interpolation 
Cubic  Interpolation 

132  mm 

417  mm 

451  mm 

445  mm 

4 -2  Measured  data  results 

The  final  results  of  this  research  were  obtained  by  applying  the  methods  dis¬ 
cussed  in  Chapter  3  to  the  downsampled  (16-by-16)  measured  data  described  in  Chap¬ 
ter  2.  The  results  of  the  methods  will  be  shown  and  compared  on  a  40-by-40  grid  in 
order  to  diminish  the  edge  noise  seen  in  hgure  2.7.  The  algorithm  ran  for  130  itera¬ 
tions  and  resulted  in  a  range  RMSE  of  721  milimeters.  The  RMSE  versus  iteration 
number  is  shown  in  hgure  4.5. 

The  interpolators  were  then  applied  to  the  downsampled  measured  data.  The 
pixel  replication  interpolator  achieved  a  range  RMSE  of  779  milimeters.  The  linear 
interpolator  resulted  in  a  range  RMSE  of  787  milimeters.  The  cubic  interpolator 
achieved  a  range  RMSE  of  788  milimeters.  The  estimated  range  graphs  of  the  inter- 
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Figure  4.3:  Range  RMSE  vs  iterations  after  applying  the  proposed  algorithm  to  the 
second  target.  This  graph  shows  a  great  reduction  in  RMSE  happening  between  zero 
and  approximately  seventy-hve.  The  algorithm  reached  approximate  convergence  at 
a  slower  rate  for  this  target  as  compared  to  the  first  target. 

polators  and  proposed  algorithm  are  shown  in  figure  4.6.  The  comparison  of  all  the 
methods  applied  to  the  measured  data  is  shown  in  table  4.3. 
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Figure  4.5:  Range  RMSE  vs  iterations  after  applying  the  proposed  algorithm  to 

the  measured  data.  The  graph  shows  the  most  RMSE  reduction  occurring  between 
zero  and  forty  iterations.  The  algorithm  reaches  convergence  after  approximately  the 
hundredth  iteration. 


Table  4.3:  Method  comparison  for  the  measured  data 


Method 

Range  RMSE 

Proposed  Algorithm 
Pixel  Replication 
Linear  Interpolation 
Cubic  Interpolation 

721  mm 

779  mm 

787  mm 

788  mm 
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Figure  4.6:  (a)  40-by-40  3-D  range  truth  of  the  measured  data.  Measured  data 

40-by-40  3-D  estimated  ranges:  (b)  Pixel  replication  estimated  range,  (c)  Linear 
estimated  range,  (d)  Cubic  estimated  range,  (e)  EM  algorithm’s  estimated  range. 


37 


V.  Conclusions  and  Future  Work 


This  section  details  conclusions  that  were  drawn  from  the  results  of  this  research. 
Future  research  potential  is  also  presented  here. 

5. 1  Conclusions 

This  research  proves  that  fusion  of  2-D  and  3-D  LADAR  images  through  EM 
increases  the  range  accuracy  of  3-D  images.  The  combination  of  2-D  high  spatial 
resolution  images  and  3-D  FLASH  LADAR  images  produces  a  new  LADAR  system 
with  improved  resolution  over  current  realizable  FLASH  3-D  sensors.  The  algorithm’s 
direct  solution  for  the  range  and  bias  allows  it  to  be  applied  to  both  measured  and 
simulated  data,  as  proved  in  this  research. 

Several  conditions  were  used  to  create  three  data  sets  for  this  research.  The 
results  for  each  data  set  show  the  EM  solution’s  range  estimation  is  a  vast  improve¬ 
ment  over  both  no-processing  and  interpolation.  This  case  is  made  clearer  with  the 
results  shown  in  the  simulated  multi-building  target  (second  target),  in  which  the 
algorithm  makes  a  sixty-five  percent  improvement  over  the  best  interpolation  results. 
The  EM  solution  was  created  under  the  assumption  that  the  2-D  and  3-D  data  was 
statistically  independent,  while  this  was  the  case  for  the  simulated  data  it  was  not 
for  the  measured  data.  In  most  cases  when  using  two  cameras  the  data  and  noise 
will  be  statistically  independent.  Given  this  finding  the  EM  algorithm  still  was  an 
improvement  over  interpolation  for  the  measured  data.  If  the  data  was  statistically 
independent  the  proposed  algorithm  would  have  done  better.  While  the  EM  solution 
may  be  more  computationally  intensive  and  require  a  second  2-D  camera,  the  range 
accuracy  would  be  a  good  trade-off  given  the  inaccuracy  of  interpolation. 

The  one  shot  capability  the  fusion  of  2-D  and  3-D  FLASH  LADAR  images 
provides  would  work  as  fast  or  faster  than  microscanning  LADAR  systems.  Again, 
LADAR  microscanning  can  involve  latency  due  to  the  fact  of  the  numerous  cubes  it 
needs.  The  algorithm  for  fusing  2-D  and  3-D  LADAR  images  opens  up  the  possibility 
of  new  LADAR  capabilities. 
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5.2  Future  Work 


The  algorithm  proved  successful  using  registered  2-D  and  3-D  images.  Future 
work  involving  a  registration  method  built  into  the  algorithm  would  further  develop 
this  type  of  system  as  a  stand  alone  system  that  includes  a  2-D  camera  and  3-D 
LADAR  system,  registration  is  not  the  only  issue  involved  with  using  two  different 
cameras.  In  most  cases  the  cameras  do  not  possess  the  same  sampling,  as  well  as 
the  images  from  both  cameras  are  uncorrelated.  Some  type  of  calibration  as  well  as 
statistical  analysis  of  the  systems  would  make  the  algorithm  work  better  and  make  a 
stand  alone  system. 

As  evidenced  in  the  measured  data,  noise  reduction  in  3-D  LADAR  systems 
would  also  make  this  algorithm  work  more  effectively  in  range  estimation.  The  mea¬ 
sured  data  was  very  noisy  and  contained  numerous  spikes  that  if  suppressed  would 
allow  the  algorithm  to  perform  better.  Estimating  noise  may  be  a  possibility  to  reduce 
noise  in  the  images.  Performing  gain  variation  compensation  and  enacting  a  photon 
counting  algorithm  would  also  improve  the  algorithm  with  respect  to  the  measured 
data. 

The  algorithm  only  considered  a  hxed  OTF,  using  blind  deconvolution  to  esti¬ 
mate  the  OTF  would  only  bolster  the  performance  of  the  algorithm.  The  work  only 
considered  range  and  bias  estimation,  the  more  estimation  built  into  it  the  better 
range  accuracy  the  system  would  be  able  to  achieve. 

Future  work  involving  a  comparison  of  the  algorithm  to  microscanning  would 
prove  what  the  best  method  of  processing  LADAR  images  would  be.  The  comparison 
should  take  into  account  the  time  it  takes  to  process  imagery  and  range  accuracy  of 
each  method.  Using  measured  data  from  LADAR  cameras  would  be  a  fairer  compar¬ 
ison  of  the  two  processing  types. 
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